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TECHNICAL MEMORANDUM 


STATISTICALLY GENERATED WEIGHTED CURVE FIT OF RESIDUAL FUNCTIONS 

FOR MODAL ANALYSIS OF STRUCTURES 

I. INTRODUCTION 


In industry, the use of analytical tools to aid in design and analysis of a product is greatly increas- 
ing. The type of tools vary depending on the type of analysis and application of the final product. In 
recent years, the method of finite element analysis (FEA), which is a mathematical representation, has 
become an industry-wide accepted tool. FEA, which is a reduction of an infinite system to a Finite one, is 
an approximation of the actual system. FEA can be applied to a variety of problems, such as stress 
analysis, thermal analysis, fluid flow, dynamics, etc. 

In the area of dynamics, an FEA can be used to determine the natural frequencies and mode 
shapes of a structural model. If the frequency of a forcing function is the same as one of the structure’s 
natural frequencies, the dynamic characteristics of the structure can produce increased elastic deforma- 
tion of the structure's material, causing fatigue, resulting in the destruction of the structure. Since the 
FEA only approximates the dynamic characteristics of the structure, it is important to eliminate any 
additional error that may be incorporated into the mathematical model. 

Errors arise in two areas: the modeling of the structure (creating the finite element model (FEM)) 
and fabricating the structure. When creating the FEM, simplifying assumptions lead to errors. In the 
fabrication of the structure, exact reproduction of the design drawings is impossible. Variations in thick- 
nesses of machined parts and limits on tolerances can change the final characteristics of the structure. 

The process of eliminating the errors in a mathematical model is known as correlation. Correlation 
involves modifying the model to match the results obtained from testing the physical structure. 

Once a mathematical model has been correlated, it can be used to predict, with increased confi- 
dence, the response of the structure under expected loadings. How a structure will behave in its opera- 
tional environment is of great importance. The structure interacts with the environment at its interface 
points. In structural dynamics, the interface points are where the object receives its excitation and 
transmits its vibrational loads back to the environment. The interface points are areas where a high level 
of confidence in the mathematical model representation is critical. 

Dynamically correlating a mathematical model with test results requires a modal test of the 
physical structure. The type of modal testing performed depends on the structure and the type of inter- 
facing the structure has with its environment. The two most commonly performed modal tests are the 
free-free and fixed-base tests. For structures with constrained (fixed) interfaces, historically the fixed- 
base modal test was performed. Fixed-base testing has many drawbacks for large structures such as 
interface simulation and cost. Extreme difficulties arise when trying to simulate the interfaces between 
the structure and its environment when only selected degrees of freedom (DOF) are to be constrained. 
The cost involved in design and construction of a test stand can become prohibitive. 

As stated previously, fixed-base tests of large structures are extremely difficult and expensive. 

As an example, the Boeing Company in Huntsville, AL, designed and built a large test stand to perform 
modal tests on the space station modules in a fixed-base configuration (fig. 1). Extreme care was taken 



ORIGINAL PAGE 

■slack and white photograph 



2 


Figure 1. Fixed-base modal test stand 





in the design of the test stand to keep the dynamic characteristics of the stand from corrupting the test 
data. A fixed-base modal test of the space station common module prototype (SSCMP), built by Boeing 
for preliminary testing, was conducted in the test stand for test stand checkout (fig. 2). The module has 
seven DOF where it interfaces with its environment. An elaborate off-loading system, which utilized air 
bags to lift the trunnions off the lower bearing surface to reduce the contact friction, was used to allow 
movement in the remaining DOF at the interface points. 

The fixed-base modal test of the SSCMP failed. The off-loading system did not free the remain- 
ing DOF at the interface locations. Boeing is presently redesigning the areas around the interface points 
of the test stand. Over a million dollars has already been spent on the test stand and its design. It can 
easily be seen why an alternative modal testing technique is greatly needed. 

Due to simulation problems and economic pressures, alternative modal testing methods are being 
considered. Mass additive and residual flexibility testing methods are two such methods. The mass 
additive method involves the attachment of masses to the interface DOF which helps to exercise the 
interface modes, aiding in constrained modes prediction (fig. 3). References 1, 2, and 3 discuss the 
implementation of the method along with its advantages and disadvantages. 

The residual flexibility modal testing method consists of measuring the free-free natural fre- 
quencies and mode shapes along with the interface frequency response functions (FRF’s). Analytically, 
the residual flexibility method has been investigated in detail, 4-6 but has not been implemented exten- 
sively for model correlation due to difficulties in data acquisition. In recent years, improvement of data 
acquisition equipment has made possible the implementation of the residual flexibility method. 7-11 The 
residual flexibility modal testing technique is applicable to a structure with distinct points (DOF) of con- 
tact with its environment. Examples are shown in figure 4. 

Payloads of the space shuttle are structures that have distinct interface points at its support trun- 
nions and are prime candidates for the residual flexibility modal test method (fig. 5). Since the 
Challenger shuttle incident in 1986, NASA requires that all final loads analysis of the payload be 
performed with a mathematical model that has been correlated to modal test data. A load analysis 
consists of coupling the payload’s mathematical model to a space shuttle mathematical model, 
dynamically exciting it analytically, and evaluating the response loads the payload exerts on the shuttle 
during takeoff, on orbit, and landing configurations. Usually the mode shapes that produce the highest 
value in a load analysis are the modes that excite the bending of the trunnions and exercise their local 
stiffness. Johnson Space Center (JSC), which is responsible for the shuttle fleet, requires the 
mathematical model to be correlated, up to 35 Hz, to a modal test data of a fixed-base configuration. A 
fixed-base configuration involves rigidly constraining the interface DOF. JSC requires this configuration 
because the exact interface point characteristics on the shuttle are unknown. These characteristics also 
change when different payloads are placed in the cargo bay. JSC also chose the fixed-base test since the 
corresponding modes are required in a loads analysis. 

If the payload was tested in a free-free configuration, the bending stiffness of the trunnions 
would not be measured (trunnion natural frequencies are usually greater than 150 Hz). Measuring the 
free-free bending stiffness of the trunnions in a free-free modal test and correlating the mathematical 
model are not sufficient to predict constrained modes. 7 In addition to measuring the bending stiffness of 
the trunnions, the residual flexibility modal testing method also measures the contributions of the higher 
order modes. 

The study of the implementation and analysis usability of the residual flexibility modal test 
technique is being conducted at the NASA/Marshall Space Flight Center (MSFC) in Huntsville, AL. 
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Figure 2. SSCMP fixed-base modal test. 



5 


Figure 3. Mass additive modal test. 
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Figure 5. Payload installed in the space shuttle. 
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A research task entitled “Residual Flexibility In-House Implementation Study” is being conducted by 
the Dynamics and Loads Branch and the Experimental Dynamics Test Branch at MSFC. During this 
study, a problem was encountered when trying to extract the residual flexibility from the test data. 

This report discusses the development of a statistically generated weighted curve fit for assisting 
in the analysis of the residual flexibility test data. First, a brief analytical discussion is presented, fol- 
lowed by an explanation of the generation of the residual function from test data, and the implementa- 
tion of the statistically generated weighted curve fit to extract the residual terms. The residual flexibility 
method is then applied to a simple beam with a trunnion type interface and a space shuttle pallet simula- 
tor. The residual terms are then compared to a mathematical model taken as the exact answer and per- 
centage errors are presented. 


II. RESIDUAL FLEXIBILITY ANALYTICAL APPROACH 


The residual technique, which was first developed by MacNeal, 5 incorporates the use of the 
higher frequency mode effects (residual modes) which are not directly used in the analysis to improve 
the reduced mathematical model. MacNeal derives statically the deflection influence coefficient matrix, 
G, associated with the free coordinates, G cc , and the redundant support coordinates, G ss , for the deter- 
minate restrained DOF, which is given by 


G = 


G cc G cs 
G Is G ss 


( 1 ) 


where G cs is the coupling of the redundant support and free coordinates. 

MacNeal further defines the deflection influence coefficient matrix with all boundary DOF 
restrained as. 


G cc = G cc~ G cs G ss G cs ■ ( 2 ) 

Next, casting the modal representation in a manner similar to the deflection influence coefficient matrix, 
equation (1), gives. 


G m = d> -K~ l Q> T 

Cl l Cl 


(3) 


where the O a term contains the eigenvectors at the connection point and Ki is the corresponding stiff- 
ness matrix. Subtracting equation (3) from equation (2) produces the residual flexibility matrix. 


CC — t'cc'Wc 


(4) 


This approach is difficult to implement the way it was derived. MacNeal only makes use of the residual 
stiffness and neglects the mass and damping terms. 

Rubin 4 starts with the same statics approximation as MacNeal, but formulates the residual flexi- 
bility in a different way. The statics approximation for a discrete system is written as. 
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(5) 


U { p=GF , 

where ) is the displacement matrix of the flexible body, F is the external force vector, and G is the 
flexibility matrix of the free body. 

The flexibility matrix is defined as the inverse of the stiffness matrix, G = K~ l . For the uncon- 
strained case, the stiffness matrix is singular and G cannot be computed directly. For the constrained 

case, G c can be calculated from K c (constrained stiffness), G c = K~ l . G and G c differ by the rigid-body 
inertia forces of the system. When the system is unconstrained, the system has elastic and rigid-body 
movement. The rigid-body movement of the system has inertia forces associated with its mass. When 
the system is constrained, the elastic motion of the system is present but the rigid-body movement is 
missing. A way to incorporate the free-body inertia into G c is needed to produce an equivalent unre- 
strained flexibility matrix. Rubin writes the displacement of the physical coordinates due to generalized 
rigid- body displacements as. 


Ur-Q>rQr , (6) 

with the differential equation relating the generalized rigid-body displacements to the external forces 
being. 


MrQr - ^rF y (7) 

T 

where M R - <& R M <I> R is the rigid-body generalized mass, and M is the physical mass matrix. The sum 
of the external forces, F, and the inertia forces, F t , can be expressed as, 

F+Fj = F-Ml) R = AF , (8) 


where 


A = . 

The constrained equation for the first-order statics approximation for a discrete system, where G c 
is the constrained-body flexibility matrix, is stated as, 

U« ) = G c (F+F i ) = G c AF . (9) 

Rubin further explains that this is an unusual statics problem where the flexibility matrix G c is singular 
and, in general, nonsymmetric. 

By adding certain generalized rigid-body displacements, Q R , to the constrained deflection given 
in equation (9), the first-order flexible-body displacements can be written as, 

uf = U»h<t> R Q R . (10) 
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Defining the in a way that it is orthogonal to the rigid-body modes (i.e., removing any rigid-body 
contributions) gives. 


4>rM[u1 1) +<I> r Q r ] = 0 . ( 11 ) 

Solving for Q R , substituting into equation (10), and using the relationship M R = oj A/ <X>/j , the free- 
body symmetric flexibility matrix, G, for the first-order approximation is, 

uf=GF , (12) 

where 

G = A t G c A , 

and A is given in equation (8). Note that G is singular since G c is singular. 

The second-order statics approximation is obtained by including the inertial and damping forces 
of the first-order deflections in equation (9), 


uf = g c a 


F-MV. ( f l) -CU ( f l) F 


(13) 


and in a similar way removing the rigid-body contributions, 

vf = GF-HF-BP 


(14) 


where 


H = GMG , 

B = GCG , 

where H and B are the inertial and damping coefficient matrices, respectively. Note that the inertial and 
damping coefficients are calculated from pre- and postmultiplication of G found from the first-order 
approximation. In addition, H and B are symmetric since G, M, and C are symmetric. 

At this point, the approximate deflection for the entire flexible body for first and second orders 
are given by equations (12) and (14), respectively. By removing the retained modes (measured modes) 
in the above equations, the residuals are obtained. 

The residuals for the first-order static equation are, 

V%=G p F , (15) 
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where 


G p =G-G N , G n = Q> N KpfQ>Jf , 

and the mode shape and generalized stiffness matrices of the retained modes are dV and K^, respec- 
tively. Similarly, the second-order formulation is expanded from the first-order statics equation by 
including correction terms. 


uf p =G p F-H p P-B p F , (16) 

where 

H p = H-H n , Hft = GftMGfl , 

B p =B-Btf , B N = GtfCGtf . 

G p is given in equation (15), H and B are given in equation (14). In addition, by making use of modal 

orthogonality, the residual mass and damping can be computed by pre- and postmultiplying the full mass 
and damping matrices, respectively. 


77p = G p M G p , 

(17) 

B p = GpCGp . 

(18) 


One now has the equations to compute the residual flexibility (stiffness), residual mass, and 
residual damping terms for a mathematical model of a structure. Note: Martinez 6 restated the equations 
to determine the residual terms in a form which is similar to the Craig-Bampton formulation. 12 The 
Craig-Bampton formulation is used widely in loads analysis. 


III. RESIDUAL FLEXIBILITY MODAL TEST APPROACH 


Modal testing of a structure is required for correlation of a mathematical model of the structure. 
A mathematical model correlated to a free-free modal test is not accurate enough to predict the con- 
strained mode shapes, as discussed in reference 7. The model must also be correlated to at least the 
residual flexibility terms, which in this study represents the flexibility around the interface points. In this 
procedure, the determinations of the test-generated Hpa, Bpa, and Gpa, the residual mass, damping, and 
flexibility, respectively, are used for comparison to the mathematical model’s residual terms found in 
equations (15) and (16). 

The residual flexibility modal test consists of first measuring the free-free modes in the fre- 
quency range of interest. This includes the rigid-body modes (in the analysis, the mathematical model’s 
rigid-body modes can be used if the total mass and its distribution in the model are relatively accurate). 
The measured free-free modes are used to correlate the global modes of the mathematical model which 
can be performed at a later time. The structure is then excited at the points where it interfaces with its 
environment (boundary or interface DOF). 
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Exciting each interface point individually, the response due to the excitation is measured at the 
same location and in the same direction as the excitation, producing a drive point FRF. Modal parame- 
ters (mode shapes and natural frequencies) are then extracted from the FPF obtained from the free-free 
modal test. From the modal parameters, synthesized FRF’s are generated and then subtracted from the 
corresponding drive point test FRF’s. The remaining functions are called the residual functions. The 
residual flexibility term is defined as the value of the residual function evaluated at zero frequency 
(constant term in the curve fit equation which will be discussed in detail later). The residual mass is the 
curvature of the residual function (second-order term in the curve fit equation) in the displacement/force 
domain. Finally, the residual damping (first-order term in the curve fit equation) is equated to the 
damping obtained from the modal test. The residual terms are then compared to the corresponding theo- 
retical residual terms. The mathematical model is modified to match the test results. Once the mathemat- 
ical model is correlated to the free-free mode shapes and the residual terms, it can be used to predict the 
constrained mode shapes. 7 

In this derivation, the residual terms will be extracted from the FRF in the displacement/force 
domain instead of the usual acceleration/force domain. Using the displacement/force domain allows 
easier recognition of the characteristics of the residual function and the residual flexibility term. 

Y a , the FRF matrix is represented in the form, 

U a {co) = Y a (co)F a (o)) . (19) 

The elements of the main diagonal of Y a {ci)) are the drive point FRF’s at each interface point and direc- 
tion. One column of Y a (co) is the result of exciting one interface point and measuring the response at 
that point and all the other interface points and directions. Y a ((d) can be thought of as a three-dimen- 
sional matrix with two dimensions of the matrix defining the location and direction of a particular FRF. 
The third dimension along the co axis represents the actual values of the FRF’s. The FRF’s are only 
measured in the frequency range of interest. The unmeasured range will produce the residual function. 

The technique for obtaining the response functions is left up to the test engineer. Care must be 
taken to accurately measure the first antiresonance of the response function. An accurate determination 
of the first antiresonance is the governing factor in determining the residual flexibility term. The first 
antiresonance is the lowest in frequency distinguishable characteristic of the FRF and, since the residual 
flexibility term is evaluated at zero frequency of the residual function, any shift in the first antiresonance 
significantly affects the residual flexibility term. 

The modal parameters are extracted from the FRF’s producing retained natural frequencies and 
mode shapes. These natural frequencies and mode shapes together with the rigid-body modes are sub- 
tracted from the corresponding FRF to produce the residual function coefficients, Y^ito ) , 

V®>= • GO) 

where <t> Ra and ® Na are the rigid-body and retained free-free modes, respectively, and. 
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y n = a~ n 1 m- n 1 , 

A N = -a) 2 +i2£ n Q)Ct)„+co 2 • 

Yr and Y n are the modal dynamic flexibility matrices for classically damped modes. The rigid-body 
mass and retained elastic-body modal mass, Mr and M^r, respectively, are identity matrices if the modes 
are initially normalized to modal mass. 

Once the residual function is determined, the residual terms (flexibility, mass, and damping) are 
extracted. A second-order form, similar to equation (16), is used to estimate the residual function, 

Y pa {(0) = (d l H pa -i(oB pa +G pa . (21) 

As can be seen, the residual function Y pa (co) is a complex number, and the real and imaginary com- 
ponents can be separated. Y pa (co) is in the form, 

Y pa (co) = a+bi , (22) 

where 

a = Real ( Y pa {co) ) , b = Imag (Y^ico) ) . 

Equating equations (21) and (22) gives, 

a+bi = (o) 2 H pa +G pa )-io)B pa . (23) 

Equating the imaginary components, the residual damping term is found to be, 

B pa = -7o ■ (24) 

The residual mass and flexibility are determined by equating the real components, 

a = G pa +co 2 H pa . (25) 

Rewriting equation (25) in matrix form gives, 

a = AX , (26) 

where 

1 a) 2 
1 co\ 

1 (D 2 




13 



and N is the total number of data points. Solving for X gives, 

X = A r'a , < 27 ) 

the residual flexibility and the residual mass. Since A is not square, other solution methods must be 
implemented, such as the pseudoinverse to determine the unknowns, G ^ and H ^ . These values are 
compared to the residual terms calculated from the mathematical model for correlation purposes. 


IV. LEAST-SQUARES CURVE FIT OF RESIDUAL FUNCTION 


The determination of the residual terms (residual mass and flexibility) from the residual function 
requires a second-order polynomial curve fit of equation (25). The imaginary term has previously been 
determined, equation (24), and is not included in the curve fit. If the subtraction of the synthesized FRF 
from the test (measured) FRF was clean (i.e., no ragged data regions), a direct curve fit could be 
accomplished for equation (27). But due to the errors in the synthesized FRF, the resulting residual 
function has regions of ragged data (fig. 6). The errors in the synthesized FRF are attributed to the limi- 
tation in the evaluation of the modal parameters from the test FRF’s. 



Figure 6. Test-generated residual function. 


A theoretical residual function in the displacement/force domain (fig. 7) has the characteristics of 
a relatively flat line in the lower frequencies and a slight upward curvature in the higher frequency 
range. These are the characteristics of a second-order equation with only the constant and second power 
terms (equation (25)). Using these characteristics, a weighted curve fit of the residual function could be 
used to eliminate the ragged data. A method of determining the weighting function for the curve fit that 
eliminates any guess work from the analyst is required. 
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Figure 7. Theoretically generated residual function. 

As stated previously, the residual function is relatively flat at the lower frequencies, i.e., the 
difference (variance) between adjacent data points is small. In the regions of ragged data, the variance 
between adjacent data is large. When the variance is large between data points, its influence on the curve 
fit should be minimized. In addition, the data point pairs with small variance should have their influence 
increased. Since a large variance requires small influence and small variance requires large influence, 
the inverse of the variance is used to eliminate the regions of ragged data. By taking the inverse of the 
variance of adjacent data points, a weighting function W, is generated. 


W{j) = 


1 n - 1 

s 2 (j) y+mi(a=l) 


x 

f) 


(x—X ) 2 


(28) 


where 


j = int . int + 1 » - ’ N ~ int ’ 

and s 2 (j) is the variance, 13 N is the total number of data points, n is the number of adjacent data points 
for which the variance is calculated, X is the mean of the n adjacent data points for each j, and int(#) is 
the integer part of # (i.e., int{ 1.5) = 1). 

Forming W into a square matrix is accomplished by setting the main diagonal equal to the 
weighting function and all other terms equal to zero. The square matrix W must be NxN in size. When 
calculating the weighting function, the first and last few data points, depending on the value of n chosen. 
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will not have a weighting value assigned. Forcing a value of zero at these points will make Wbe NxN. 
This can be justified due to the fact that at the extreme limits, data acquisition equipment is not as accu- 
rate as in the middle range of the equipment, and the test generated residual function is not as accurate 
there either. Eliminating these points will have little influence on the curve fit. This is demonstrated in 
the study on the beam in the following sections with exaggerated data point elimination. In addition, 
only the weighting values of about two or four data points will be set equal to zero compared to about 
250 total data points. 

Premultiplying both sides of equation (26) by the weighting function matrix gives, 

Wa = WAX . (29) 

Multiplying both sides by A 1 and solving for X yields, 

X = (A T WAr l A T Wa . (30) 

Since A T WA is a square matrix, the inverse can be computed directly. In this form, there is no need to 
normalize the weighting function. The residual terms (residual mass and flexibility), which are contained 
in X, can be used to correlate the mathematical model. 


V. APPLICATION TO BEAM WITH APPENDAGE 

One of the simplest structures to model is a straight beam. A theoretical solution for the natural 
frequencies and mode shapes of an unrestrained (free-free) uniform beam is available. 14 The natural fre- 
quencies are given by, 


fi = 



/= 1,2,3, ... 


(31) 


where 


^ 1.2. 3.4.5 = 4.73004074, 7.85320462, 10.9956078, 14.1371655, 17.2787597 , 


A i = (2/+l)|; i > 6, 


and the mode shapes are given by, 


"(f) 


cosh 


X : x 


+ cos 


X;X 


cr, sinh 


X:X 


+ sin 


X ; x 


( 32 ) 


where 


^ 1 . 2 . 3 . 4.5 = 0.982502215, 1.000777312,0.999966450, 1.000001450,0.999999937 , 

Cj = 1.0 ; / >6 . 
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With the “exact” answer available, the mathematical model’s elastic modes can be correlated. Once the 
model is correlated, the model is said to be exact. Table 1 shows the frequency difference less than one- 
half percent and modal assurance criteria (MAC) values of 1.00. The MAC is an averaged value for the 
comparison of the mode shapes in the range from zero to one. A value of one indicates a high degree of 
similarity between mode shapes. 

Table 1. Correlation of beam model and theoretical solution. 


Mode 

Number 

Test 

Frequency (Hz) 

Model 

Frequency (Hz) 

Difference 

(Hz) 

Difference 

(Percent) 

MAC 

1 

9.9764 

9.9740 

0.024 

0.02 

1.00 

2 

27.4715 

27.493 

-0.0215 

0.07 

1.00 

3 

53.8007 

53.899 

-0.0983 

0.18 

1.00 

4 

88.8484 

89.097 

-0.2486 

0.27 

1.00 

5 

132.6023 | 

133.096 

-0.4937 

0.37 

1.00 

6 

185.0520 

185.895 

-0.8448 

0.45 

1.00 


The correlated model can be used as a baseline for methodology, program development, and test 
method evaluation. 

For this study, the beam structure was modified by attaching a trunnion simulator to one end (fig. 
8). The trunnion simulator was a brass rod with two aluminum plates attached on both ends. The plates 
were used to connect the brass rod to the beam and to attach the testing instruments. The trunnion simu- 
lator gives the beam a distinct stiffness change from the free-free beam. The mass of the trunnion simu- 
lator is small relative to the mass of the beam and does not affect correlation to the solution of the theo- 
retical uniform free-free beam. 

The free-free modal test hardware setup consists of a Hewlett Packard 3562 for data acquisition. 
The accelerometer utilized during the model test was a PCB Piezotronics 303A accelerometer. The sys- 
tem was excited with an impact hammer utilizing a PCB Piezotronics 208A03 load cell. Once the 
measurements were acquired, they were ported to a Hewlett Packard 9000/750 workstation running 
Leuven Measurement Systems (LMS) MODAL™ package for modal parameter estimation. For ease of 
testing, the accelerometer was fixed at the interface point and the excitation hammer was moved. 

Initially, the beam was tested in a free-free configuration (fig. 9), and the elastic modes were 
compared to the mathematical model. Table 2 shows the percentage difference between natural fre- 
quencies and the MAC values. The measurements and the modal parameter estimation (natural frequen- 
cies and mode shapes) of the free-free modal test are verified by the results in table 2. 

By exciting the test structure at the point where it interfaces with its environment (the end of the 
brass rod) and measuring the response at the same location, a drive point FRF in the displacement 
domain was obtained (fig. 10). The modal parameters were computed for the first six natural frequen- 
cies, and a synthesized response function was generated (fig. 1 1). The first six natural frequencies in the 
maximum frequency range of 0 to 200 Hz were chosen arbitrarily. 
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Figure 8. Beam with trunnion attachment test article. 
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Figure 9. Free-free modal test of beam. 


Table 2. Comparison of beam model and test data. 


Mode 

Number 

Test 

Frequency (Hz) 

Model 

Frequency (Hz) 

Difference 

(Hz) 

Difference 

(Percent) 

MAC 

1 

9.9888 

9.9764 

0.0124 

0.12 

1.00 

2 

27.4918 

27.4715 

0.0203 

0.07 

1.00 

3 

53.8265 

53.8007 

0.0258 

0.04 

1.00 

4 

88.9338 

88.8484 

0.0854 

0.09 

1.00 

5 

132.5651 

132.6023 

-0.0372 

0.02 

1.00 

6 

184.7864 

185.0520 

-0.2656 

0.14 

1.00 


For visual comparison, the drive point and synthesized FRF’s have been plotted together in 
figure 12. The frequency range has been extended to show the relationship between the measured fre- 
quencies and the first bending frequency of the trunnion, which occurs at 289 Hz. In an actual test, 
everything above the cutoff frequency would be unmeasured. The unmeasured region of the FRF is 
responsible for the residual curve. The subtraction of the test and synthesized response functions produc- 
ing the residual function can be seen in figure 13. The characteristics of the theoretical residual function 
are apparent, but with regions of ragged data produced by modal parameter estimations. 

A direct curve fit of the residual function was made (fig. 14). The curve fit did not converge on 
the underlying second-order characteristic curve. Figure 14 was plotted with a linear vertical axis since 
the value of the curve becomes negative around 150 Hz and the characteristics of the second-order fit 
can be seen. The ragged data are orders of magnitude higher than the underlying second-order character- 
istic curve, which throws off the curve fit. By stepping through the data points, a weighting value 
(inverse of variance) with respect to the neighboring data points can be calculated (equation (28)). The 
neighboring data points were initially examined in sets of three. It can be seen in figure 15, which con- 
tains the weighting and residual functions, that the weighting function has low values where the residual 
function has regions of ragged data. Another curve fit was performed using the statistically generated 
weighting function (equation (30)) with highly accurate results (fig. 16). 

The exact answers were obtained by equation ( 16). The same values can be obtained by follow- 
ing the test procedure in generating the residual function. The drive-point FRF is generated by using the 
elastic modes ot the model. Figure 17 shows the residual function generated from the mathematical 
model. The residual values of 1 ,6045e-5 for flexibility and 5.7540e-12 for mass were compared to the 
exact answer (mathematical model) and found to have 0.81- and 18.44-percent error, respectively. It has 
been determined from the research task being performed by MSFC/NASA that the residual mass and 
damping terms can be neglected in the correlation of the mathematical model, which is used to predict 
constrained natural frequencies and mode shapes. However, the residual mass term is important in the 
curve fitting process due to the curvature of the residual function. To verify the process, additional 
weighting functions and curve fits were performed by varying parameters. The number of adjacent data 
points was changed to two and four along with the range of the data points for which the weighting 
functions were computed and applied (table 3). 

By varying the upper and lower frequency ranges, an evaluation of the effects of excluding data 
points in the curve fit process is examined. The changes in the percent error of the frequency range 
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Figure 10. Test drive-point FRF for beam 
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Figure 1 1. Synthesized drive-point FRF for beam 
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Figure 12. Overlay of test and synthesized drive-point FRF’s for beam 
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Figure 13. Residual function. 
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Frcquency(Hz) 

Figure 14. Direct curve fit of residual function. 
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Figure 15. Weighting function for beam 
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Frequency(Hz) 

Figure 16. Weighted curve fit of residual function for beam. 
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Figure 17. Theoretically generated residual function from beam model. 




Table 3. Residual flexibility curve fitting errors for beam. 


Lower 

Frequency 

Upper 

Frequency 

2 Points 

Percent 

Error 

3 Points 

Percent 

Error 

4 Points 

Percent 

Error 

Exact 


1.6165e-5 

0.00 

1.6165e-5 

0.00 

1.6165e-5 

0.00 

0 

200 

1.5761e-5 

2.52 

1.6035e-5 

0.81 

1.5833e-5 

2.06 

30 

200 

1.5757e-5 

2.53 

1.6035e-5 

0.80 

1.5835e-5 

2.04 

60 

200 

1.577 le-5 

2.44 

1.6225e-5 

0.37 

1.5865e-5 

1.86 

0 

180 

1.5758e-5 

2.52 

1.6066e-5 

0.61 

1.5888e-5 

1.71 

30 

180 

1.5758e-5 

2.52 

1.6067e-5 

0.60 

1.5890e-5 

1.70 

60 

180 

1.5772e-5 

2.43 

1.6268e-5 

0.64 

1.5923e-5 

1.50 


variation in table 3 is insignificant. The comparison of the utilization of two, three, and four consecutive 
data points to generate the weighting function shows a small change with three points producing the best 
results. 


The larger error in the two- and four-point evaluation can be attributed to the distribution of the 
data point of the curve. For instance, if two consecutive data points straddle a peak (define the peak) at 
equal heights, their variance may be low but their overall magnitude is greater than the characteristic 
curve. A high weighting value would be assigned that would produce an incorrect curve fit (i.e., residual 
flexibility value). 

After the test mode shapes and natural frequencies have been obtained, the data manipulation 
was performed by programs written in MATLAB™ by The MathWorks. The program MK-Phi-FRF.m 
takes the mathematical model’s mass and stiffness matrices and generates the natural frequencies, mode 
shapes, and appropriate drive-point FRF for the response function approach. MK-Phi-FRF.m sets up the 
input parameters for the program rbtcfv.m. The program rbfcfv.m, which utilizes the response function 
approach, computes the residual function from the drive-point FRF and generated elastic and rigid-body 
response tunctions. A weighting function is generated and applied to the curve fit of the residual func- 
tion, producing the residual terms. The program res.m uses the matrix approach to generate the residual 
terms from a mathematical model. The programs MK-Ph-FRF.m, rbfcfv.m, and res.m can be found in 
the appendices. Each program is documented and has the appropriate equation reference. 


VI. APPLICATION TO PALLET SIMULATOR 


Once the method and programs have been verified on a simple structure (beam) with a known 
exact (theoretical) answer, the implementation on a more complex structure is to be performed. A frame 
structure (fig. 18) has been designed and built to simulate a space shuttle cargo pallet. The fundamental 
mode distribution and trunnion modes of the frame structure are similar to a space shuttle cargo pallet. 

A theoretical “exact” solution is not available for the pallet simulator. To obtain the closest 
model representation to the physical structure, the mathematical model (FEM) is to be correlated to the 
elastic modes obtained from the modal test data. All the elastic modes up to the first bending mode of 
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Figure 18. Space shuttle pallet simulator. 





the trunnions are to be correlated. In an actual test, only a subset of the elastic modes (up to a prechosen 
cutoff frequency) would be used in the correlation. In this analysis, the first eight natural frequencies are 
used in the determination of the residual values. Since the mathematical model is correlated through the 
14 elastic modes, the use of the model to predict the residual terms of the first eight modes can be 
accomplished with confidence. 

The free-free modal test hardware setup consists of a Hewlett Packard 9000/750 workstation 
with a Hewlett Packard 3565S frontend for data acquisition (fig. 19). The accelerometers utilized during 
the modal test were PCB Piezotronics STRUCTCEL™ accelerometers. LMS software was employed, 
with their FMON™ package for data acquisition and the MODAL™ package for modal parameter 
estimation. 

The pallet simulator was suspended by bungy cords in a free-free environment. The modal test 
was conducted using random vibration shaker excitation in the three translational directions in sequence 
(fig. 20). The mathematical model was then correlated to the first 14 elastic natural frequencies and 
mode shapes. The results of the correlation are found in table 4. With a frequency difference less than 
1 percent and MAC values of 1.00, an “exact” mathematical model is said to be obtained. 

The drive-point FRF’s are to be measured at the points where the structure would connect 
(interface) with its environment (end of the trunnions). A hammer impact excitation was used to excite 
the structure with a data acquisition frequency range up to 100 Hz. If the shaker was attached to the end 
of the trunnion, mass and stiffness loading would occur due to the size of the trunnion and would impose 
additional uncertainties in the data. The response function at the other interface points was acquired for 
each drive-point excitation with a total of 144 FRF’s obtained. 

Examining one of the drive-point FRF’s, modal parameters (mode shapes and natural frequen- 
cies) were extracted for the first eight elastic modes. Using these eight modes, a synthesized FRF was 
generated. The synthesized and the drive-point FRF’s are plotted together (fig. 21) for visual compari- 
son. The frequency range in figure 21 has been extended to show their relationship. Subtracting the 
synthesized FRF from the drive-point FRF, the residual function was generated (fig. 22). It can be seen 
that a direct curve fit of the residual function would cause extreme divergence as in the beam’s residual 
function curve fit. 

A weighting function was generated using three consecutive data points (fig. 23). Applying the 
weighted curve fit to the residual function, the residual flexibility and residual mass are found to be 
7.6738e-6 and 5.6104-12, respectively (fig. 24). The residual terms from the “exact” mathematical 
model were computed and compared to the test-generated residual terms with test errors of 2.50- and 
1.36-percent for flexibility and mass, respectively. The number of consecutive data points and the 
frequency range for weighting function generation were changed to verify the curve fitting process 
(table 5). Similar variation in the percent error was observed as in the beam results. 
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Figure 19. Modal test setup and test equipment. 
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Figure 20. Free-free modal test configuration. 
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Frequency (Hz) 

Figure 22. Residual function for pallet simulator. 
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Figure 23. Weighting function for pallet simulator. 
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Figure 24. Weighted curve fit of residual function for pallet simulator. 





Table 4. Correlation results for pallet simulator. 


Mode 

Number 

Test 

Frequency (Hz) 

Model 

Frequency (Hz) 

Difference 

(Hz) 

Difference 

(Percent) 

MAC 

1 

16.5620 

16.5274 

-0.0346 

0.21 

1.00 

2 

21.2079 

21.3203 

0.1124 

0.53 

1.00 

3 

46.4941 

46.1600 

-0.3341 

0.72 

1.00 

4 

51.3995 

51.5989 

0.1994 

0.39 

1.00 

5 

75.7266 

75.3772 

-0.3494 

0.46 

1.00 

6 

83.7987 

84.4266 

0.6279 

0.75 

1.00 

7 

96.7800 

96.1199 

-0.6601 

0.68 

1.00 

8 

98.6837 

98.4312 

-0.2525 

0.26 

1.00 

9 

105.7185 

106.0403 

0.3218 

0.30 

1.00 

10 

112.8742 

113.1375 

0.2633 

0.23 

1.00 

11 

135.6785 

136.2627 

0.5842 

0.43 

1.00 

12 

136.6981 

137.9833 

1.2852 

0.94 

1.00 

13 

138.6954 

139.4763 

0.7809 

0.56 

1.00 

14 

167.3594 

167.7916 . 

0.4322 

0.26 

1.00 


Table 5. Curve fitting errors of pallet simulator residual flexibility. 


Lower 

Frequency 

Upper 

Frequency 

2 Points 

Percent 

Error 

3 Points 

Percent 

Error 

4 Points 

Percent 

Error 

Exact 


7.8707e-6 


7.8707e-6 


7.8707e-6 


0 

100 

7.5776e-6 

3.72 

7.6738e-6 

2.50 

7.6923e-6 

2.27 

20 

100 

7.6056e-6 

3.37 

7.6842e-6 

2.37 

7.7037e-6 

2.12 

40 

100 

7.8214e-6 

0.63 

7.9305e-6 

0.76 

7.9457e-6 

0.95 

0 

90 

7.2536e-6 

7.84 

7.5488e-6 

4.09 

7.5553e-6 

4.01 

20 

90 

7.3175e-6 

7.03 

7.5658e-6 

3.87 

7.5764e-6 

3.74 

40 

90 

7.7639e-6 

1.36 

7.9176e-6 

0.60 

7.9368e-6 

0.84 
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VII. CONCLUSION 


An applicable least-squares curve fitting procedure using a statistically generated weighting 
function has been developed. This curve fitting procedure enables the accurate curve fitting of test- 
generated residual functions that contain regions of ragged data. The ragged data is the result of the 
inability to accurately evaluate the modal parameters from the FRF during the modal testing process. 
Programs for calculating the residual terms from modal test data and a mathematical model (FEM) have 
been written in MATLAB™. The curve fit of the drive-point residual function for a beam with a 
trunnion type interface, retaining the first six elastic modes, produced residual flexibility values with 
errors in the 0.4- to 2.5-percent range. The errors were determined by comparing the residual terms of 
the test data to the residual terms generated from a mathematical model which was correlated to the 
exact theoretical solution. 

The curve fitting process was also applied to a frame structure that simulates a space shuttle 
cargo pallet, which is more complex than a straight beam. The residual flexibility values for the pallet 
simulator, retaining the first eight elastic modes, were found to be in the 0.6- to 4.0-percent range. A 
theoretical solution is not available for the frame structure. The “exact” answer was generated by corre- 
lating the mathematical model to all the elastic natural frequencies and mode shapes up to the first 
bending mode of the trunnion (i.e., 14 elastic modes). 

The development of the weighted curve fit of the residual function has made possible further 
investigation of the residual flexibility modal testing technique. Additional studies are needed in the area 
of modal parameter estimation, since this is the cause of the ragged data in the residual function. The 
next step, which is also under investigation, is how to use the residual flexibility to accurately correlate a 
mathematical model. 
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APPENDIX A 


PROGRAM MK-Phi-FRF.m 


The program MK_Phi_FRF.m takes the mathematical model’s mass and stiffness matrices and 
generates the natural frequencies, modes shapes, and the appropriate drive-point frequency response 
function for the response function approach. MK_Phi_FRF.m sets up the input parameters for the pro- 
gram rbfcfv.m for the mathematical model data. The input parameters are set up to run the beam prob- 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


PROGRAM MK * Phi * FRF . m 

™ iS £™ grain COmputs from the and Stiffness matrices 

of a FEM program, the eigenvalues and eigenvectors. It then 
computes the Frequency Resopnce Function for a specified 

^° n , Ce P ° lnt and citation point. The frequencies and 
mode shapes are splat up into rigid body and elastic 
components . 


Note: The mass & stiff 


matrices must be stored in matlab 


———'WWW i(IU ^ L JJT 

format with the matrices names as "M r 
respectivily . 


"K' 


clear 
echo off 

prename=input ( 'Enter PreName? 


I Se )j f °f foiiowing generating FRF corresponding to 

% ?" Fu11 ' 1 -Elastic, 2-Rigid Body, 3-No FRF 

Flag=0; 


3=26; 
k=26; 
nrbm= 6 ; 
nrem=6 ; 
xdof s= [27] ; 
StFreq=0 . 0; 
step=0 . 78125 ; 
dim=512 ; 


% drive coordinate 
% responce coordinate 
% number of rigid body modes 
% number of elastic modes 
% DOF to exclude from analysis 
% start freqency in Hz. 

% delta-f, freqency step in Hz. 
% size of frf to be generated 


eval ( [ ' load ' , prename, ' .MK.mat ' ] ) ; 
[m,m] =size (M) ; 

FRF=zeros ( 1 , dim) ' ; 
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f=linspace(StFreq,dim*step,dim) ' ; 

f =f (1 : dim) ; 

Rad=f * (2*pi ) ; 

% compute mode shape and frequencies 

disp(' Computing Eigenvalues and Eigenvectors') 

[ Phi , w, Freq, GM] =EignSol2 (M, K) ; 

% separate elastic and rigid body modes 

disp(' Separate Elastic and Rigid Body Modes') 

Phi (xdof s , :)=[]; 

Phie=Phi ( : ,nrbm+l :nrbm+nrem) ; 

Phirb=Phi ( : , l:nrbm) ? 

Freqe=Freq (nrbm+1 : nrbm+nrem) ; 

Freqrb=Freq ( 1 : nrbm) ; 

% save modal parameters 

eval ( [ ’ save ' , [prename , ' . Elas . mat J , Freqe Phie ] ) 
eval ( [ ' save ', [prename, ' .RB. mat' ], ' Freqrb Phirb']) 

% Generate Full FRF 
if Flag==0 

disp{' Proccessing Full FRF') 

Lambda=zeros (dim, m-1 ) ; 

Phi2=Phi ( j , : ) • *Phi (k, : ) ; 

Rad2=Rad. ~2; 
r=l; 

while r<=m-l 

FRF=FRF+ ( Phi 2 (r) . / (Freq(r) ~2-Rad2) ) ; 
r=r+l ; 
end 
end 

% Generate Elastic FRF 
if Flag==l 

disp ( ' Proccessing Elastic FRF') 

Lambda=zeros (dim, nrem) ; 

Phi2=Phie ( j , : ) .*Phie(k, : ) ; 

i=l; 
r=l ; 

while r<=nrem 
while i<=dim; 

Lambda (i , r) = (Phi2 (r) / (Freqe (r) ~2-Rad ( i ) '2 ) ) ; 

i=i+l ; 
end 
i=l ; 

FRF ( : ) =FRF ( : ) +Lambda ( : , r ) ; 
r=r+l ; 
end 
end 

% Generate Rigid Body FRF 
if Flag==2 



disp(' Proccessing Rigid Body FRF') 
xf = zeros (dim, nrbm) ; 

Phi 2= Phi rb ( j , : ) . *Phirb(k, : ) ; 

i = l; 
r=l ; 

while r<=nrbm 
while i<=dim; 

Lambda (i , r ) = { Phi 2 ( r ) / (Freqrb ( r ) "2 -Rad ( i ) "2 ) ) ; 
i = i + l ; 
end 
i = l ; 

FRF ( : ) =FRF ( : ) +Lambda ( : , r ) ; 
r=r+l ; 
end 
end 


% save FRF in (displacement/f ) format 

eval ( [ 'save ' , [prename, ' . ' ,int2str(j) , ' . ' ,int2str (k) , 
' .FRF. mat' ] , ' FRF' ] ) ; 
disp(' FRF Saved in x/f format') 


% plot FRF 

semi logy ( f , abs (FRF) ) 

title ([' Filename ', prename, ' : Drive num2str (j ) , ... 

' - Response num2str (k) ] ) 
xlabel ( ' Freq' ) 
ylabel ( 'X/F' ) 

% DISCRI PTION OF TERMS 

drive point DOF 
response point DOF 
number of rigid body modes 
number of elastic modes 
DOFS to exclude from the analysis 
frequency for the FRF to start 
frequency step (resolution) in Hz. 
dimension of the FRF (size) 
vector containing the FRF 
frequency values in Hz. 
frequency values in radians 
mass matrix 
stiffness matrix 
eigenvectors 
eigenvalues in radians 

eigenvalues in radians on main diagonal 
generalized mass 

elastic and rigid body eigenvectors 


% 1 
% k 

% nrbm 
% nrem 
% xdofs 
% StFreq 
% step 
% dim 
% FRF 
f 

% rad 
% M 
% K 
% Phi 
% Freq 
% w 
% GM 

% Phie,Phirb 


% 


% Freqe, Freqrb elastic and rigid body eigenvalues 


43 




APPENDIX B 


PROGRAM rbf.m 


The program rbfcfv.m, which utilizes the response function approach, computes the residual 
function from the drive-point FRF and generated elastic and rigid-body response functions. A weighting 
function is generated and applied to the curve fit of the residual function producing the residual terms. 
The input parameters are set up to run the beam problem. 


% PROGRAM rbf.m 
% 

% Computes the Residual Function and performs a weighted least 
% square curve fit of the Residual Function, 
clear 
echo off 
i=sqrt (-1) ; 


% 


Required input 
3 = 26 ; 
k=2 6 ; 

HzStep=0. 78125; 
n=512 ; 

HzStFreq=0. 78125; 
dptl=100 ; 
dptu=256 ; 


num ' sample=3 ; 


% drive point DOF 
% response point DOF 
% squared resolution (Delta Hz) 

% number of data points 
% start frequency in Hertz 
% lower data pt to consider in curve fit 
% upper data pt to consider in curve fit 
% dptl Sc dptu for flexibility of range 
% in curve fit 

% # of data samples to use in forming weighting 
% function (range 2-6). Note: if 5 or 6 change 
% w(r+l,r+l)= to w(r+2,r+2)= 


% Load Data from MATLAB Format Files 

% prename . Elas .mat contains Free-Free Frequencies and Modes 
% ( "Freqe, Phie" ) 

% prename . FRF .mat contains Test FRF {"FRF") 

% prename . RB . mat contains Rigid Body Frequencies and Modes 
% ("Freqrb, Phirb" ) 

prename=input { ' Enter PreName? ' , ' s ' ) ; 
eval ( [ ' load ' , prename, ' .Elas. mat' ] ) ; 
eval { [ 'load ' , prename, ' . FRF. mat ' ] ) ; 
eval ( [ 'load ' , prename, ' . RB.mat ' ] ) ; 


PttSC&Diti'J rVb,A i. AUH. HO { f IvJVXS 
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% Setup Parameters and Initialize Matrices 
[tmp,nrem]=size(Phie) ; 

[ tmp # nrbm] =size (Phirb) ; 

Hz =1 inspace (HzStFreq, n*HzStep, n) ' ; 

Hz=Hz (1 :n) ; 

RadSqd= (Hz . * ( 2*pi ) ) . ~2; 

FRF=FRF ( 2 : n+1 ) ; 

FRFe=zeros ( 1 , n) # ; 

FRFrb=zeros (1 , n) ' ; 
w=zeros (n, n) ; 
xx=zeros (2,1) ; 

% Synthesize Elastic Mode FRF, from Equ . 3.2 
disp(' Synthesize Elastic Mode FRF') 

Phie2=Phie ( j , : ) . *Phie (k, : ) ; 
r=l ; 

while r<=nrem 

FRFe=FRFe+ (Phie2 (r) . / (Freqe (r) "2 . -RadSqd) ) ; 
r=r+l ; 
end 

% Generate Rigid Body FRF, from Equ. 3.2 
disp(' Generate Rigid Body FRF ' ) 

Phirb2=Phirb ( j , : ) . *Phirb (k, : ) ; 
r=l ; 

while r<=nrbm 

FRFrb=FRFrb+ (Phirb2 (r ) . / (Freqrb(r) "2 . -RadSqd) ) ? 
r=r+l ; 
end 

% Subtract Ridged Body and Elastic FRFs from Test FRF, 

% from Equ . 3.2 

retFRF= (FRFrb) + (FRFe) ; 

ResFRF= ( (FRF) + (retFRF) ) ; 

% Computing Weighting Function, Equ. 4.1 
disp(' Computing Weighting Function') 
b=abs (ResFRF ) ; 
r=dptl ; 

while r<=dptu- (num * sample-1 ) 
w(r+l , r+1) =1/ ( std(b(r : r+ (num ‘ sample- 1 ) ) ) ) "2 ; 
if w(r+l,r+l) == Inf 

disp ( ' *** Variance = 0.0 Substituting 1.0e+10 ***') 
w (r+1 , r+1) =1 . 0e+*10; %value may need to be changed to fit mean 
end 
r = r + l ; 
end 

% Normalize Weighting Function 
w= w . /max ( max ( w ) ) ; 

% Forming "A" matrix, from Equ. 3.8 
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A=RadSqd; 

A ( : , 2 ) =ones (n, 1) ; 


% Solving for Unknowns, from Equ. 4.3 
disp(' Solving for Unknowns') 
xx= (A' *w*A) \ (A' *w*b) ; 

Gr=xx (2,1) ; 

Hr=xx (1,1) ; 

% Computing Curve Fit, from Equ. 3.8 
disp{' Computing Curve Fit') 
cfit=A*xx; 

% Saving Residual Terms and FRF 

eval ( [ 'save ', [prename Res .mat '], ' Gr Hr ResFRF']); 

% Plot Residual FRFs 

semi logy (Hz f abs (ResFRF) , ' r- ' , Hz , abs ( retFRF) , ' g- ' , . . . 

Hz, abs (FRF) , 'w: ' ) 

title (' Frequency Response Functions for Residual Flexibility') 
xlabel ( 'Frequency (Hz) ' ) 
ylabel ( 'Displacement (x/f ) ' ) 
legend ( 'Residual Function' , 'Synthesized' , 'Measured' ) ; 

%print 

disp(' ** Screen is PAUSED. Enter or Return to continue ** ' ) 
pause 

% Plot Curve Fit ant Etc. 

semi logy (Hz , abs (b) , Hz , abs (cfit ) , Hz , abs (diag (w) ) ) 
title {' Difference FRF & Curve Fit') 
xlabel ( 'Frequency (Hz) ' ) 
ylabel ( 'Displacement (x/f) ' ) 

legend ( 'Residual Function' , 'Curve Fit ', 'Weighting Function'); 


ResTerms= . . . 

sprintf('\n Residual Terms\n Gr = %10.8e\n Hr = %10 . 8e ' , Gr , Hr ) 
disp (ResTerms ) ; 


% DESCRIPTION OF TERMS 


% A 
% cfit 
% FRF 
% FRFe 
% FRFrb 
% Gr , Hr 
% Hz 
% nrem 
% nrbm 

% Phie,Freqe 
% Phie2 

% Phirb,Freqrb 
% Phird2 
% RadSqd 


matrix of coefficients of the unknowns 
curve fit of residual function 
test frequency response function 
synthesized FRF from elastic modes 
synthesized FRF from rigid coefficient modes 
residual flexibility and residual mass 
vector of frequencies in hertz for plotting 
number of retained elastic modes 
number of rigid body modes 
elastic modal properties 
the square of each term in Phie 
rigid body modal properties 
the square of each term in Phirb 
vector of frequencies in radians, terms squared 



% ResFRF 
% retFRF 
% w 
% xx 


residual function 

synthesized FRF of retained modal parameters 
matrix with weighting function on main diagonal 
contains residual mass and flexibility 
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APPENDIX C 


PROGRAM res.m 


The program res.m uses Rubin's matrix approach to generate the residual terms from a mathe- 
matical model. This program also computes the constrained natural frequencies and mode shapes from a 
subset of free-free modes and drive-point residual flexibilities . 7 The input parameters are set up to run 

thP hpcim nmhlpm * 


% PROGRAM res.m 
% 

% This program computes the residuals and the constrained mode 
% shapes generated from the residual flexibility. 

% See discription of terms at end of program, 
clear 
echo off 
% input data 

prename=input ( 'Enter PreName? ','s')- 
iDof = [26 ] ; 

ConDof = [26] ; 

rm= [12345678]; 

nrbm=2 ; 

eval ( [ 'load ' , prename, ' .MK.mat' ] ) ; 

% reorder M&K so interface DOF's are in last Rows & Columns 
% generate mode Shapes ' 

disp(' Calculating eigenvalues and eigenvectors') 
niDof=length (iDof ) ; 
nConDof= length (ConDof) ; 
nrm=length(rm) ; 

[n, n] =size (M) ; 
ol= [1 ;n] ; 
ol (iDof) = [ ] ; 
ol (n-niDof+1 :n) =iDof ; 

M=M (ol , ol ) ; 

K=K (ol , ol ) ; 

[Phi,w,Freq,GM] =EignSol2 (M,K) ; 

% construct transf omation from constrained to free-free 
% assumes rigid body modes are in first columns, from Equ. 2.8 
Phirb=Phi ( : , 1 :nrbm) ; 

A=eye (n) -M*Phirb*Phirb' ; 
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% generate constrained stiffness, constrained flexibility, and 
% free-free flexibility matrix, from Equ. 2.12 
disp(' Generating flexiblity matrices') 

Kc=K; 

StConDof =n-nConDof +1 ; 

Kc ( S tConDo f : n , S tConDo f : n ) =Kc ( S tConDo i : n , S tConDo f : n ) + . . . 

diag (diag (Kc ( StConDof : n, StConDof :n) ) ) ; 

Gc=inv(Kc) ; 

G=A' *Gc*A; 


% generate retained flexibility matrix and 
% remove from full flexibility matrix, from Equ. 2.15 
disp(' Generate Retained Flexibility Matrix') 

Kn=w(rm, rm) . ~ 2 ; 

Kn=Kn (nrbm+1 : nrm, nrbm+1 : nrm) ; 

Phin=Phi ( : , rm) ; 

Phin=Phin ( : , nrbm+1 : nrm) ; 

Gn= (Phin/Kn) *Phin' ; 

Gr=G-Gn; 

% partition mode shapes and retained flexibility 
disp ( ' Partition Matrices ' ) 

Phiki=Phi ( 1 :n-niDof , rm) ; 

Phikb=Phi (n-niDof +1 :n, rm) ; 

Grib=Gr { 1 : n-niDof , n-niDof +1 : n) ; 

Grbb=Gr (n-niDof + 1 :n,n-niDof + l :n) ; 

% generate full residual mass and extract boundary mass terms, 
% from Equ . 2.17 
Hr=Gr ' *M*Gr; 

Hrbb=Hr (n-niDof + 1 :n, n-niDof + 1 :n) ; 

% THE FOLLOWING IS AN IMPLEMENTATION OF THE PROCEDURE FOUND 
% IN REFERENCE [7] . 

% generate transformation matrix 

disp(' Generate Transformation Matrix') 
tmp=Grib/Grbb; 

[ tmpm, tmpn] =size ( tmp) ; 

T=Phiki-tmp*Phikb; 

[tm, tn] =size (T ) ; 

T ( 1 : tm, tn+1 : tn+tmpn) =tmp; 

T ( tm+1 : n , 1 : tn) =zeros (n-tm, tn) ; 

T (tm+1 :n, tn+1 : tn+tmpn) =eye (n-tmpm, tmpn) ; 

% form constrained model, solve, and convert solution to 
% physical coordinates 

disp(' Form Constrained Model, Solve', and') 
disp(' Convert to Physical Coord.') 

Mb=T ' *M*T ; 

Kb=T ' *K*T ; 

Mnn=Mb ( 1 : nrm , 1 : nrm ) ; 

Knn=Kb ( 1 :nrm, 1 :nrm) ; 
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[Phinn, wc] =EignSol (Mnn, Knn) ; 

Phinn (nrm+1 :nrnw-tmpn, 1 : nrm) =zeros ( tmpn, nrm) ; 

Phic=T*Phinn; 

% exact solution by FEM 

disp(' Exact Solution by FEM') 

Ke=K ( 1 : n-nConDof , 1 : n-nConDof ) ; 

Me=M ( 1 : n-nConDof , 1 : n-nConDof ) ; 

[Phie, we, Freqe,GMe] =EignSol2 (Me,Ke) ; 

Phie {n-nConDof +1 :n, 1 :n-nConDof) = zeros (nConDof, n-nConDof ) ; 

% checking orthoginality 
Orthog= Phie ' *M* Phie ; 

% saving data 

eval ( [ ' save ’ , [prename, ' . res .mat ' ] , ' Grbb Hrbb Orthog . . 
/ascii ' ] ) ; 

disp ( ' That's All Folks') 


% DISCRIPTION OF TERMS 

% niDof number of interface DOFs . 

% iDof row or column corresponding to the interface DOF . 

% nConDof number of constrained DOF for determinatly 

% constrained model . 

% ConDof constrained DOF for determinatly constrained model 

% nrm number of retained modes. 

% rjn mode number of retained modes including r-b mode. 

% nrbm number rigid body modes . 

% m,K mass and stiffness of full model 

^ reorder to place interface DOF's in last col & rows 

% phi full mode shapes 

% w frequency in matrix form of full model 


% Phirb 
% A 
% 

% Kc 
% Gc 
% G 

% Kn , Gn 
% Gr 

% Phiki, Phi kb 
% Grib, Grbb 
% T 

% Mb, Kb 
% Mnn , Knn 
% 

% Phinn 
% wc 
% Phie 
% Me , Ke 
% 

% Phie 
% we 

% Orthog 


rigid body mode shapes 

transformation martix for constrained flexibility 
matrix to free-free flexibility matrix 
constrained stiffness matrix 
constrained flexibility matrix 
free-free flexibility matrix 
retained stiffness and flex, matrices 
residual flexibility matrix 

retained modes, interior and boundary DOF's 
subset of Gr, interior and boundary DOF's 
transformation to reduce structures coordinates 
reduced mass and stiffness matrices 
subsets of Mb and Kb respect, obtained by setting 
boundary displacements equal to zero, 
constrained mode shapes in reduced coordinates 
freq. in matrix form of constrained model 
constrained mode shapes in full coordinates 
mass and stiffness matrices for exact solution of 
constrained model 
exact constrained mode shapes 

exact freq. in matrix form for constrained model 
orthogonality of exact and constrained mode shapes 
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